/*

This file makes Appendix C Table 2

*/
global writeout 0 // change to 1 for writing out, 0 to not
/* global quietly quietly // comment out for noisy */


$quietly { // top top quietly
clear
set more off
set seed 1337
set sortseed 1337

capture cd "~/Desktop/RESTAT_CODE"

global dir_root = "."
global dir_code = "${dir_root}/code"
global dir_data = "${dir_root}/data"
global dir_output "${dir_root}/output"

// Load the bootstrap data
quietly do "${dir_code}/1.1_cache-load_lasso_data.do"
quietly do "${dir_code}/2.4_set_bootstrap_se_scalars.do"

local prefix_1 "boostrap_results_"
local prefix_2 "boostrap_results_knock_"
local prefix_3 "boostrap_results_letters_"

quietly ///
forvalues cpi = 1/3 {
  foreach cp in "c" "p" {
    use "${dir_output}/bootstrap/`prefix_`cpi''`cp'-1.dta", clear
    noi di "`cp'_`cpi'"
    foreach v in "lambda" "u_a" "u_b" "alpha_1" "alpha_0" {
      replace T_star`v' = abs(T_star`v')
      _pctile T_star`v', p(0.0000000001 95)

      scalar `v'_ci5_`cp'_`cpi' = r(r2) - r(r1)
      scalar `v'_ci95_`cp'_`cpi' = r(r2) - r(r1)
      noisily display "    `v': [`=`v'_ci5_`cp'_`cpi'' <--> `=`v'_ci95_`cp'_`cpi'']"
    }
  }
}

use "${dir_data}/table2_lasso_results.dta", clear

******************************************************************************
 * Create statistics and table
******************************************************************************

  global if1 ""
  global if2 "& (phase1_inperson == 1)"
  global if3 "& (phase2 == 1)"

  svyset geoid_blkgroup [pweight = w_together_smooth], strata(stratum_phase1) || house_id

  noisily mean cia_pos [pw = w_together_smooth ]
  noisily mean pcr_pos [pw = w_together_smooth ] if (sampled_individual ==1)
  noisily mean pcr_pos [pw = w_together_smooth ] if (sampled_individual ==1) & phase1_inperson == 1
  noisily mean pcr_pos [pw = w_together_smooth ] if (sampled_individual ==1) & phase2 ==1
  noisily mean any_pos [pw = w_together_smooth ]

  mat table_matrix = J(`=5 * 4', 3,.) // rows: 5 vars * 4 rows each (est, se, ci5, ci95), cols: 3
    mat colnames table_matrix="All Observations" "In Person" "Letters"
    mat rownames table_matrix="a_1-est" "a_1-se" "a_1-ci5" "a_1-ci95" ///
                               "a_0-est" "a_0-se" "a_0-ci5" "a_0-ci95" ///
                               "lambda-est" "lambda-se" "lambda-ci5" "lambda-ci95" ///
                               "mu_a-est" "mu_a-se" "mu_a-ci5" "mu_a-ci95" ///
                               "mu_b-est" "mu_b-se" "mu_b-ci5" "mu_b-ci95"

  quietly ///
  forvalues col = 1/3{
    local weight "[iw = w_together_smooth]"
    local symp est_smp_p_`col'
    local symc est_smp_c_`col'

    preserve
      keep if (sampled_individual ==1) ${if`col'}

      sum `symc' `weight' if cia_pos == 1
        scalar alpha_1=r(mean)
        mat table_matrix[1, `col'] = round(`=alpha_1' ,.0001) // alpha_1 estimate
        mat table_matrix[2, `col'] = round(`=(alpha_1_ci5_c_`col' + alpha_1_ci95_c_`col')/2/1.96' ,.0001) // alpha_1 SE
        mat table_matrix[3, `col'] = round(`=alpha_1 - alpha_1_ci5_c_`col'' ,.0001) // alpha_1 CI 5%
        mat table_matrix[4, `col'] = round(`=alpha_1 + alpha_1_ci95_c_`col'' ,.0001) // alpha_1 CI 95%

      sum `symc' `weight' if cia_pos == 0
        scalar alpha_0=r(mean)
        mat table_matrix[5, `col'] = round(`=alpha_0' ,.0001) // alpha_0 estimate
        mat table_matrix[6, `col'] = round(`=(alpha_0_ci5_c_`col' + alpha_0_ci95_c_`col')/2/1.96' ,.0001) // alpha_0 SE
        mat table_matrix[7, `col'] = round(`=alpha_0 - alpha_0_ci5_c_`col'' ,.0001) // alpha_0 CI 5%
        mat table_matrix[8, `col'] = round(`=alpha_0 + alpha_0_ci95_c_`col'' ,.0001) // alpha_0 CI 95%

      scalar lambda = `=alpha_1 / alpha_0'
        mat table_matrix[ 9, `col'] = round(`=lambda' ,.0001) // lambda estimate
        mat table_matrix[10, `col'] = round(`=(lambda_ci5_c_`col' + lambda_ci95_c_`col')/2/1.96' ,.0001) // lambda SE
        mat table_matrix[11, `col'] = round(`=lambda - lambda_ci5_c_`col'' ,.0001) // lambda CI 5%
        mat table_matrix[12, `col'] = round(`=lambda + lambda_ci95_c_`col'' ,.0001) // lambda CI 95%

      *odds ratio -- mu
      sureg (cia_pos `symp', nocons) (tested_yn `symp', nocons) (p_it `symp', nocons)
        cap nlcom [cia_pos]`symp'*[tested_yn]`symp' / [p_it]`symp'
        mat temp3 = r(b)
        scalar stemp3 = temp3[1,1]
        mat table_matrix[13, `col'] = round(`=stemp3' ,.0001) // tab B1 mu estimate
        mat table_matrix[14, `col'] = round(`=(u_b_ci5_p_`col' + u_b_ci95_p_`col')/2/1.96' ,.0001) // mu SE
        mat table_matrix[15, `col'] = round(`=stemp3 - u_b_ci5_p_`col'' ,.0001) // mu CI 5%
        mat table_matrix[16, `col'] = round(`=stemp3 + u_b_ci95_p_`col'' ,.0001) // mu CI 95%

    restore
  }
  mat list table_matrix

scalar row0="$\alpha^1$"
scalar row1="$\alpha^0$"
scalar row2="$\lambda$"
scalar row3="$\mu$"

quietly foreach row of numlist 0/3 {
  // Row 1 - estimate
  noi di "`=row`row'' ", _c
  foreach col of numlist 1/3 {
    noi di " & " %4.3f table_matrix[`row'*4 + 1,`col'], _c
  }
  noi di "\\"

  // Row 2 - Standard Error
  foreach col of numlist 1/3 {
    noi di "& (" %4.3f table_matrix[`row'*4 + 2,`col'] ")", _c
  }
  noi di "\\"

  // Row 3 - CI
  foreach col of numlist 1/3 {
    noi di "& [" %4.3f table_matrix[`row'*4 + 3,`col'] "," %4.3f table_matrix[`row'*4 + 4,`col'] "]", _c
  }
  noi di "\\ \addlinespace"
}

} // end top top quietly
